% Solves dynamic stochastic general equilibrium using Klein's methodadjcost
% Money growth shock
% dynsolveklein.m calls dyneqklein.m, jacob_reiter.m, and kleinsolv.m.

global Params par_;
format long
jacstep = par_.jacstep;                   % size of step for numerical jacobian calculation
eqcutoff = 10*jacstep;                    % stable-unstable eigenvalue cutoff 

phiMAT   = [phiR 0; 0 phiA];              % Matrix of shock persistence parameters 

% X = [Pdist, DP, mlag, V, C, PI, m]
nsj = 3;                                  % number of scalar jump variables (C, PI, m)
nss = 2;                                  % number of scalar state variables (DP, mlag)
nz  = 2;                                  % number of exogenous shock processes (money, TFP)
nPhi = nump;                              % number of points of the distribution (states)
nV   = nump;                              % number of points of the value function (jumps)
Ntot = nPhi+nss+nV+nsj;                   % total number of variables            

% STORE VARIABLES IN PARAMS
Params.nsj = nsj;
Params.nss = nss;
Params.nz = nz;
Params.nPhi = nPhi;
Params.nV = nV; 

% Y WILL CONSIST OF FOUR PARTS
Params.ix     = 1:Ntot;                   % length Ntot
Params.ixnext = Ntot+1:2*Ntot;            % length Ntot
Params.iz     = 2*Ntot+1:2*Ntot+nz;       % length nz
Params.iznext = 2*Ntot+nz+1:2*Ntot+2*nz;  % length nz

Params.beta = beta; 
Params.gam = gam; 
Params.chi = chi; 
Params.theta= theta; 
Params.mu = mu; 
Params.nu = nu;
Params.eta = eta;
Params.omega = omega;

Params.Rss = Rss; 
Params.Cbar = Cbar; 
Params.phiPI = NaN; 
Params.phiC = NaN; 

Params.model = model; 
Params.menucost = menucost;          
Params.lbar = lbar; 
%Params.ksi = ksi;              
Params.noise = noise;
Params.alpha = alpha;
Params.alphapos = alphapos;

Params.PMAT = PMAT; 
Params.T2 = T2; 
Params.Pgrid = Pgrid; 
Params.pstep = pstep; 

Xss = [Pdist(:);PD;mbar;Vf(:);Cbar;mu;mbar]; 
stst2=[Xss;Xss;zeros(nz,1);zeros(nz,1)];

resid = dyneqklein(stst2);
if max(abs(resid))>100*sqrt(eps)
  disp('WARNING: Large residual at stst: max(abs(resid)) = ');
  disp(max(abs(resid)));
  keyboard
end
tic

fprintf('\n\n')  
disp('START JACOBIAN CALCULATION')
   jac = jacobian(@dyneqklein,resid,stst2,jacstep);
   ajac = jac(:,Params.ixnext);
   bjac = jac(:,Params.ix);
   cjac = jac(:,Params.iznext);
   djac = jac(:,Params.iz);

fprintf('\n\n')  
disp('START KLEIN SOLUTION')

[JUMPS,STATEDYNAMICS,stableeigs,unstableeigs] = kleinsolve(ajac,bjac,cjac,djac,phiMAT,nPhi+nss,eqcutoff);
clear jac ajac bjac cjac djac